#######################
## Replication File ###
## Smoke and Mirrors ##
## May 2025 ###########
#######################

setwd("Data")

## Read in public and elite survey data
P <- read.csv("Public_Aug2024.csv")
E <- read.csv("Elite_Nov2024.csv")

## Ensure control is the held out condition
P$treatment <- factor(P$treatment, levels = c("control", "deny", "validity", "applicatory"))
E$treatment <- factor(E$treatment, levels = c("control", "deny", "validity", "applicatory"))

## function to cluster SE
cl <- function(fm, cluster){ 
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  if(length(cluster) != length(fm$fitted.values)){cluster <- cluster[-fm$na.action]}
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }


#############
## Figures ##
#############

## Figure 2: Effect of Mitigation Strategies on Certainty, Outrage

P_noIO <- P[P$IOtreatment==0,]

Certainty <- lm(certainty ~ treatment, data=P_noIO)
CS <- cl(Certainty, P_noIO$ResponseId)

certainty_data <- data.frame(
  Estimate = coef(Certainty)[-1], SE = CS[-1, 2],
  CI_low = coef(Certainty)[-1] - 1.96 * CS[-1, 2], 
  CI_high = coef(Certainty)[-1] + 1.96 * CS[-1, 2],
  Category = factor(c("Denial", "Validity Challenge", "Applicability Challenge"),
                    levels = c("Denial", "Validity Challenge", "Applicability Challenge")),
  Outcome = "Certainty"
)

Outrage <- lm(outrage ~ treatment, data=P_noIO)
OS <- cl(Outrage, P_noIO$ResponseId)

outrage_data <- data.frame(
  Estimate = coef(Outrage)[-1], SE = OS[-1, 2],
  CI_low = coef(Outrage)[-1] - 1.96 * OS[-1, 2], 
  CI_high = coef(Outrage)[-1] + 1.96 * OS[-1, 2],
  Category = factor(c("Denial", "Validity Challenge", "Applicability Challenge"),
                    levels = c("Denial", "Validity Challenge", "Applicability Challenge")),
  Outcome = "Outrage"
)

plot_data_combined <- rbind(certainty_data, outrage_data)

library(ggplot2)
library(dplyr)

combined_plot <- ggplot(plot_data_combined, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_rect(aes(
    ymin = as.numeric(as.factor(Category)) - 0.5,
    ymax = as.numeric(as.factor(Category)) + 0.5,
    xmin = -Inf, xmax = Inf, fill = Category), alpha = 0.1) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, linewidth = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = NULL, x = "", y = NULL) +
  facet_grid(rows = vars(Category), cols = vars(Outcome), scales = "free_y", space = "free") +
  scale_x_continuous(
    limits = c(-0.5, 0.42), 
    breaks = round(seq(-0.4, 0.4, 0.2), 2)) +
  scale_y_continuous(
    breaks = 1:length(levels(plot_data_combined$Category)),  
    labels = c("Denial", "Validity \nChallenge", "Applicability \nChallenge")) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    strip.text.x = element_text(size = 14, face = "bold", margin = margin(b = 15)),  
    axis.text.y = element_text(size = 14, hjust = 0),  
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank(),  
  )

combined_plot


## Figure 3: Effect of Mitigation Strategies on Willingness to Punish

Punish <- lm(punish_PC ~ treatment, data=P_noIO)
PS <- cl(Punish, P_noIO$ResponseId)

punish_data <- data.frame(
  Estimate = coef(Punish)[-1], SE = PS[-1, 2],
  CI_low = coef(Punish)[-1] - 1.96 * PS[-1, 2],
  CI_high = coef(Punish)[-1] + 1.96 * PS[-1, 2],
  Category = factor(c("Denial", "Validity Challenge", "Applicability Challenge"),
                    levels = c("Denial", "Validity Challenge", "Applicability Challenge"))
)

punish_plot <- ggplot(punish_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_rect(aes(
    ymin = as.numeric(as.factor(Category)) - 0.5,
    ymax = as.numeric(as.factor(Category)) + 0.5,
    xmin = -Inf, xmax = Inf, fill = Category), alpha = 0.1) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Willingness to Punish", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +  # Only rows are faceted
  scale_x_continuous(
    limits = c(-0.41, 0.2), 
    breaks = round(seq(-0.4, 0.42, 0.2), 2)) +
  scale_y_continuous(
    breaks = 1:length(levels(punish_data$Category)),  
    labels = c("Denial", "Validity \nChallenge", "Applicability \nChallenge")) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 14, hjust = 0),  
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

punish_plot


## Figure 4: Effect of UN intervention on willingness to punish

P$treatment_expand <- relevel(factor(P$treatment_expand), ref="control")
Punish_expand <- lm(punish_PC ~ treatment_expand, data=P)
PES <- cl(Punish_expand, P$ResponseId)

# Denial (left panel)
deny_data <- data.frame(
  Estimate = coef(Punish_expand)[7:9], SE = PES[7:9, 2],
  CI_low = coef(Punish_expand)[7:9] - 1.96 * PES[7:9, 2],
  CI_high = coef(Punish_expand)[7:9] + 1.96 * PES[7:9, 2],
  Category = factor(c("Denial", "Denial + UN", "Denial + UN (neg prime)"))
)

deny_plot <- ggplot(deny_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Denial, UN Response \n on Willingness to Punish", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +
  scale_x_continuous(
    limits = c(-0.41, 0.205), 
    breaks = round(seq(from=-0.4, to=0.2, by=0.1),2)) +
  scale_y_continuous(
    breaks = seq(1, length(levels(deny_data$Category)), by = 1),  # Adjust breaks based on compressed spacing
    labels = c("Denial", "Denial + UN", "Denial + UN\n(Neg Prime)"),
    expand = expansion(mult = c(0.02, 0.02))  # Compress space between categories
  ) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

deny_plot


# Validity (middle panel)
validity_data <- data.frame(
  Estimate = coef(Punish_expand)[10:12], SE = PES[10:12, 2],
  CI_low = coef(Punish_expand)[10:12] - 1.96 * PES[10:12, 2],
  CI_high = coef(Punish_expand)[10:12] + 1.96 * PES[10:12, 2],
  Category = factor(c("Validity", "Validity + UN", "Validity + UN (neg prime)"))
)

validity_plot <- ggplot(validity_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Validity Challenge,\n UN Response on Willingness to Punish", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +
  scale_x_continuous(
    limits = c(-0.41, 0.205), 
    breaks = round(seq(from=-0.4, to=0.2, by=0.1),2)) +
  scale_y_continuous(
    breaks = seq(1, length(levels(validity_data$Category)), by = 1),  # Adjust breaks based on compressed spacing
    labels = c("Validity\nChallenge", "Validity\nChallenge + UN", "Validity\n Challenge +\n UN (Neg Prime)"),
    expand = expansion(mult = c(0.02, 0.02))  # Compress space between categories
  ) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

validity_plot

# Applicability (right panel)
applicatory_data <- data.frame(
  Estimate = coef(Punish_expand)[2:4], SE = PES[2:4, 2],
  CI_low = coef(Punish_expand)[2:4] - 1.96 * PES[2:4, 2],
  CI_high = coef(Punish_expand)[2:4] + 1.96 * PES[2:4, 2],
  Category = factor(c("Applicability", "Applicability + UN", "Applicability + UN (neg prime)"))
)

applicatory_plot <- ggplot(applicatory_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Applicability Challenge, \n UN Response on Willingness to Punish", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +
  scale_x_continuous(
    limits = c(-0.47, 0.205), 
    breaks = round(seq(from=-0.4, to=0.2, by=0.1),2)) +
  scale_y_continuous(
    breaks = seq(1, length(levels(applicatory_data$Category)), by = 1),  # Adjust breaks based on compressed spacing
    labels = c("Applicability\nChallenge", "Applicability\nChallenge + UN", "Applicability\nChallenge + UN \n(Neg Prime)"),
    expand = expansion(mult = c(0.02, 0.02))  # Compress space between categories
  ) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

applicatory_plot


## Figure 5: Treatment Effects among Foreign Policy Elites

E$treatment_expand <- relevel(as.factor(E$treatment_expand), ref="control")

## All Elites (left panel)
Treats <- lm(punish_PC ~ treatment_expand, data=E)
coefs <- coef(Treats)[-1]
ses <- cl(Treats, E$ResponseId)[2:7,2]

plot_data <- data.frame(
  Estimate = coefs, SE = ses,
  CI_low = coefs - 1.96 * ses,
  CI_high = coefs + 1.96 * ses,
  Category = c(
    "Applicability", "Applicability + UN message",
    "Denial", "Denial + UN message",
    "Validity", "Validity + UN message"),
  Panel = c(
    "Applicability Challenge", "Applicability Challenge",
    "Denial", "Denial",
    "Validity Challenge", "Validity Challenge")
)

plot_data$Category <- factor(
  plot_data$Category,
  levels = c(
    "Denial + UN message", "Denial",
    "Validity + UN message", "Validity",
    "Applicability + UN message", "Applicability"),
  labels = c(
    "Denial +\nUN message", "Denial",
    "Validity \nChallenge +\nUN message", "Validity \nChallenge",
    "Applicability \nChallenge +\nUN message", "Applicability \nChallenge")
)

plot_data$Panel <- factor(
  plot_data$Panel,
  levels = c("Denial", "Validity Challenge", "Applicability Challenge")
)

main_efx_plot <- ggplot() +
  geom_rect(
    data = plot_data,
    aes(
      ymin = as.numeric(as.factor(Category)) - 0.5,  
      ymax = as.numeric(as.factor(Category)) + 0.5,  
      xmin = -Inf,
      xmax = Inf,
      fill = Panel
    ),
    alpha = 0.1
  ) +
  # Add point ranges
  geom_pointrange(
    data = plot_data,
    aes(
      y = as.numeric(as.factor(Category)),  # Original unscaled positions
      x = Estimate,
      xmin = CI_low,
      xmax = CI_high
    ),
    shape = 16, size = 0.8
  ) +
  geom_vline(xintercept = 0, color = "#780000", size = 0.75) +
  facet_grid(rows = vars(Panel), scales = "free_y", space = "free") +
  labs(
    title = "All Foreign Policy Elites",
    x = "Effect on Willingness to Punish",
    y = NULL
  ) +
  scale_x_continuous(
    limits = c(-1.05, 0.67)  
  ) +
  scale_y_continuous(
    breaks = 1:length(levels(plot_data$Category)),  
    labels = levels(plot_data$Category)  
  ) +
  theme_minimal() +
  theme(
    strip.text = element_blank(),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
    axis.text.y = element_text(size = 14, hjust = 0, margin = margin(r = 10)),  
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90"),
    axis.line = element_blank(),
    legend.position = "none"
  )

main_efx_plot


## Global North elites (middle panel)

EGN <- E[E$CountryType=="Global North",]
EGS <- E[E$CountryType=="Global South",]

Treats <- lm(punish_PC ~ treatment_expand, data=EGN)
coefs <- coef(Treats)[-1]
ses <- cl(Treats, EGN$ResponseId)[2:7,2]

plot_data <- data.frame(
  Estimate = coefs,
  SE = ses,
  CI_low = coefs - 1.96 * ses,
  CI_high = coefs + 1.96 * ses,
  Category = c(
    "Applicability", "Applicability + UN message",
    "Denial", "Denial + UN message",
    "Validity", "Validity + UN message"),
  Panel = c(
    "Applicability Challenge", "Applicability Challenge",
    "Denial", "Denial",
    "Validity Challenge", "Validity Challenge")
)

plot_data$Category <- factor(
  plot_data$Category,
  levels = c(
    "Denial + UN message",
    "Denial",
    "Validity + UN message",
    "Validity",
    "Applicability + UN message",
    "Applicability"
  ),
  labels = c(
    "Denial +\nUN message",
    "Denial",
    "Validity \nChallenge +\nUN message",
    "Validity \nChallenge",
    "Applicability \nChallenge +\nUN message",
    "Applicability \nChallenge"
  )
)

plot_data$Panel <- factor(
  plot_data$Panel,
  levels = c("Denial", "Validity Challenge", "Applicability Challenge")
)

GN_efx_plot <- ggplot() +
  geom_rect(
    data = plot_data,
    aes(
      ymin = as.numeric(as.factor(Category)) - 0.5,  
      ymax = as.numeric(as.factor(Category)) + 0.5,  
      xmin = -Inf,
      xmax = Inf,
      fill = Panel
    ),
    alpha = 0.1
  ) +
  geom_pointrange(
    data = plot_data,
    aes(
      y = as.numeric(as.factor(Category)),  
      x = Estimate,
      xmin = CI_low,
      xmax = CI_high
    ),
    shape = 16, size = 0.8
  ) +
  geom_vline(xintercept = 0, color = "#780000", size = 0.75) +
  facet_grid(rows = vars(Panel), scales = "free_y", space = "free") +
  labs(
    title = "Global North Elites",
    x = "Effect on Willingness to Punish",
    y = NULL
  ) +
  scale_x_continuous(
    limits = c(-1.05, 0.67)  
  ) +
  scale_y_continuous(
    breaks = 1:length(levels(plot_data$Category)),  
    labels = levels(plot_data$Category) 
  ) +
  theme_minimal() +
  theme(
    strip.text = element_blank(),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),  
    axis.text.y = element_text(size = 14, hjust = 0, margin = margin(r = 10)),  
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90"),
    axis.line = element_blank(),
    legend.position = "none"
  )

GN_efx_plot

## Global South Elites (right panel)
Treats <- lm(punish_PC ~ treatment_expand, data=EGS)
coefs <- coef(Treats)[-1]
ses <- cl(Treats, EGS$ResponseId)[2:7,2]

plot_data <- data.frame(
  Estimate = coefs,
  SE = ses,
  CI_low = coefs - 1.96 * ses,
  CI_high = coefs + 1.96 * ses,
  Category = c(
    "Applicability", "Applicability + UN message",
    "Denial", "Denial + UN message",
    "Validity", "Validity + UN message"),
  Panel = c(
    "Applicability Challenge", "Applicability Challenge",
    "Denial", "Denial",
    "Validity Challenge", "Validity Challenge")
)

plot_data$Category <- factor(
  plot_data$Category,
  levels = c(
    "Denial + UN message",
    "Denial",
    "Validity + UN message",
    "Validity",
    "Applicability + UN message",
    "Applicability"
  ),
  labels = c(
    "Denial +\nUN message",
    "Denial",
    "Validity \nChallenge +\nUN message",
    "Validity \nChallenge",
    "Applicability \nChallenge +\nUN message",
    "Applicability \nChallenge"
  )
)


plot_data$Panel <- factor(
  plot_data$Panel,
  levels = c("Denial", "Validity Challenge", "Applicability Challenge")
)

GS_efx_plot <- ggplot() +
  geom_rect(
    data = plot_data,
    aes(
      ymin = as.numeric(as.factor(Category)) - 0.5,
      ymax = as.numeric(as.factor(Category)) + 0.5,
      xmin = -Inf,
      xmax = Inf,
      fill = Panel
    ),
    alpha = 0.1
  ) +
  geom_pointrange(
    data = plot_data,
    aes(
      y = as.numeric(as.factor(Category)), 
      x = Estimate,
      xmin = CI_low,
      xmax = CI_high
    ),
    shape = 16, size = 0.8
  ) +
  geom_vline(xintercept = 0, color = "#780000", size = 0.75) +
  facet_grid(rows = vars(Panel), scales = "free_y", space = "free") +
  labs(
    title = "Global South Elites",
    x = "Effect on Willingness to Punish",
    y = NULL
  ) +
  scale_x_continuous(
    limits = c(-1.05, 0.67), 
  ) +
  scale_y_continuous(
    breaks = 1:length(levels(plot_data$Category)),
    labels = levels(plot_data$Category)  
  ) +
  theme_minimal() +
  theme(
    strip.text = element_blank(),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),  
    axis.text.y = element_text(size = 14, hjust = 0, margin = margin(r = 10)),  
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90"),
    axis.line = element_blank(),
    legend.position = "none"
  )

GS_efx_plot



#################################
## Appendix Tables and Figures ##
#################################

##################
## Appendix A.1 ##
##################

## Table A4: Public Sample Demographics

round(prop.table(table(P$pid)), 2)
round(prop.table(table(P$age_cat)),2)
round(prop.table(table(P$edu_cat)),2)
round(prop.table(table(P$race)),2)
round(prop.table(table(P$hisp)),2)
round(prop.table(table(P$income)),2)
round(prop.table(table(P$region_cat)),2)


## Table A5: Effect of Mitigation Strategies on Certainty, Outrage

Certainty <- lm(certainty ~ treatment_expand, data=P)
Outrage <- lm(outrage ~ treatment_expand, data=P)

cl_cert <- cl(Certainty, P$ResponseId)
cl_outrage <- cl(Outrage, P$ResponseId)

library(stargazer)
stargazer(Certainty, Outrage,
          se = list(cl_cert[,"Std. Error"], cl_outrage[,"Std. Error"]),
          dep.var.labels = c("Certainty", "Outrage"),
          covariate.labels = c("Applicability Challenge", "Denial", "Validity Challenge"),
          omit = c("UN", "Constant"),
          omit.stat = c("f", "ser", "rsq")
)

## Table A6: Effect of Mitigation Strategies on punishment

WTP <- lm(punish_PC ~ treatment_expand + regimetype + pastbehavior, data=P)
cl_WTP <- cl(WTP, P$ResponseId)

stargazer(WTP,
          se = list(cl_WTP[,"Std. Error"]),
          dep.var.labels = "Willingness to Punish",
          covariate.labels = c("Applicability Challenge", "Denial", "Validity Challenge",
                               "Democratic State", "History of Compliance"),
          omit = c("UN", "Constant"),
          omit.stat = c("f", "ser", "rsq")
)


## Table A7: Effect of Mitigation Strategies + UN responses on punishment

Punish_expand <- lm(punish_PC ~ treatment_expand, data=P)
PES <- cl(Punish_expand, P$ResponseId)

stargazer(Punish_expand, 
          se=list(PES[,"Std. Error"]),
          dep.var.labels = "Willingness to Punish",
          covariate.labels = c("Applicability Challenge", "Applicability Challenge + UN Response",
                               "Applicability Challenge + UN Response, neg prime", "Denial",
                               "Denial + UN Response", "Denial + UN Response, neg prime",
                               "Validity Challenge", "Validity Challenge + UN Response",
                               "Validity Challenge + UN Response, neg prime"),
          omit = c("Constant", "control"),
          omit.stat = c("f", "ser", "rsq"))


## Table A8: Government of Employment in Elite Survey Sample

employers <- sort(prop.table(table(E$GovEmployer[E$GovEmployer != ""])), decreasing=T)
round(employers, 2)[1:25]


## Table A9: Treatment Effects on Punishment, Foreign Elites

Treats.All <- lm(punish_PC ~ treatment_expand, data=E)
ses.All <- cl(Treats.All, E$ResponseId)

Treats.GN <- lm(punish_PC ~ treatment_expand, data=EGN)
ses.GN <- cl(Treats.GN, EGN$ResponseId)

Treats.GS <- lm(punish_PC ~ treatment_expand, data=EGS)
ses.GS <- cl(Treats.GS, EGS$ResponseId)

stargazer(Treats.All, Treats.GN, Treats.GS,
          se=list(ses.All[,2], ses.GN[,2], ses.GS[,2]),
          dep.var.labels = "Willingness to Punish",
          covariate.labels = c("Applicability Challenge", "Applicability Challenge + UN Response",
                               "Denial", "Denial + UN Response", "Validity Challenge", 
                               "Validity Challenge + UN Response"),
          omit = c("Constant"),
          omit.stat = c("f", "ser", "rsq"))


## Table A10: Willingness to Punish (Avg) and Willingness to Impose Sanctions

WTP_avg <- lm(punish_mean ~ treatment_expand + regimetype + pastbehavior, data=P)
avg.SES <- cl(WTP_avg, P$ResponseId)

WTP_san <- lm(punish_econ ~ treatment_expand + regimetype + pastbehavior, data=P)
avg.econ <- cl(WTP_san, P$ResponseId)

stargazer(WTP_avg, WTP_san,
          se=list(avg.SES[,2], avg.econ[,2]),
          dep.var.labels = c("Avg. Punishment", "Sanctions"),
          covariate.labels = c("Applicability Challenge", "Applicability Challenge + UN Response",
                               "Applicability Challenge + UN Response, neg prime", "Denial",
                               "Denial + UN Response", "Denial + UN Response, neg prime",
                               "Validity Challenge", "Validity Challenge + UN Response",
                               "Validity Challenge + UN Response, neg prime", 
                               "Democrtic Transgressor", "History of Compliance"),
          omit = c("Constant", "control"),
          omit.stat = c("f", "ser", "rsq"))


## Table A11: Willingness to Punish by issue area

WTPA <- lm(punish_mean ~ treatment_expand, data=P[P$issue=="aggression",])
WTPA.SES <- cl(WTPA, P$ResponseId[P$issue=="aggression"])

WTPO <- lm(punish_mean ~ treatment_expand, data=P[P$issue=="human rights",])
WTPO.SES <- cl(WTPO, P$ResponseId[P$issue=="human rights"])

stargazer(WTPA, WTPO,
          se=list(WTPA.SES[,2], WTPO.SES[,2]),
          dep.var.labels = "Willingness to Punish",
          covariate.labels = c("Applicability Challenge", "Applicability Challenge + UN Response",
                               "Denial", "Denial + UN Response", "Validity Challenge", 
                               "Validity Challenge + UN Response"),
          omit = c("Constant", "control", "neg prime"),
          omit.stat = c("f", "ser", "rsq"))


## Table A12: Het FX of UN Intervention

P_UN <- P[grepl("control", P$treatment_expand) | grepl("+ UN", P$treatment_expand),]

WTP <- lm(punish_PC ~ treatment_expand, data=P_UN)
cl_WTP <- cl(WTP, P_UN$ResponseId)

WTP_NP <- lm(punish_PC ~ treatment_expand, data=P_UN[P_UN$IOprime=="negative",])
names(WTP_NP$coefficients) <- gsub(", neg prime", "", names(WTP_NP$coefficients))
cl_WTP_NP <- cl(WTP_NP, P_UN$ResponseId[P_UN$IOprime=="negative"])

WTP_LC <- lm(punish_PC ~ treatment_expand, data=P_UN[P_UN$UN_confidence_score==0,])
cl_WTP_LC <- cl(WTP_LC, P_UN$ResponseId[P_UN$UN_confidence_score==0])

WTP_HC <- lm(punish_PC ~ treatment_expand, data=P_UN[P_UN$UN_confidence_score==1,])
cl_WTP_HC <- cl(WTP_HC, P_UN$ResponseId[P_UN$UN_confidence_score==1])

P_UN$Resp_Norms_Score <- ifelse(P_UN$Resp_Norms >=4, 1, 0)

WTP_LN <- lm(punish_PC ~ treatment_expand, data=P_UN[P_UN$Resp_Norms_Score==0,])
cl_WTP_LN <- cl(WTP_LN, P_UN$ResponseId[P_UN$Resp_Norms_Score==0])

WTP_HN <- lm(punish_PC ~ treatment_expand, data=P_UN[P_UN$Resp_Norms_Score==1,])
cl_WTP_HN <- cl(WTP_HN, P_UN$ResponseId[P_UN$Resp_Norms_Score==1])

stargazer(WTP, WTP_NP, WTP_LC, WTP_HC, WTP_LN, WTP_HN,
          se=list(cl_WTP[,2], cl_WTP_NP[,2], cl_WTP_LC[,2], 
                  cl_WTP_HC[,2], cl_WTP_LN[,2], cl_WTP_HN[,2]),
          dep.var.labels = "Willingness to Punish",
          covariate.labels = c("Applicability Challenge + UN",
                               "Denial + UN", 
                               "Validity Challenge + UN"),
          #keep = c("deny +", "applicatory + "), 
          omit = c("neg", "control", "Constant"),
          omit.stat = c("f", "ser", "rsq"))


## Table A13: How certainty & outrage predict punishment, Global South vs. Global North

CO <- lm(punish_PC ~ certainty + outrage, data=E)
se.CO <- cl(CO, E$ResponseId)

COGN <- lm(punish_PC ~ certainty + outrage, data=EGN)
se.COGN <- cl(COGN, EGN$ResponseId)

COGS <- lm(punish_PC ~ certainty + outrage, data=EGS)
se.COGS <- cl(COGS, EGS$ResponseId)

CO2 <- lm(punish_PC ~ certainty*CountryType + outrage*CountryType, data=E)
se.CO2 <- cl(CO2, E$ResponseId)

stargazer(CO, COGN, COGS,
          se=list(se.CO[,2], se.COGN[,2], se.COGS[,2]),
          dep.var.labels = "Willingness to Punish",
          omit = c("Constant"),
          omit.stat = c("f", "ser", "rsq"))



## Figure A1: Effect of Mitigation Strategies on Punishment (Avg Measure)

Punish <- lm(punish_mean ~ treatment, data=P_noIO)
PS <- cl(Punish, P_noIO$ResponseId)

punish_data <- data.frame(
  Estimate = coef(Punish)[-1], SE = PS[-1, 2],
  CI_low = coef(Punish)[-1] - 1.96 * PS[-1, 2],
  CI_high = coef(Punish)[-1] + 1.96 * PS[-1, 2],
  Category = factor(c("Denial", "Validity Challenge", "Applicability Challenge"),
                    levels = c("Denial", "Validity Challenge", "Applicability Challenge"))
)

punish_plot <- ggplot(punish_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_rect(aes(
    ymin = as.numeric(as.factor(Category)) - 0.5,
    ymax = as.numeric(as.factor(Category)) + 0.5,
    xmin = -Inf, xmax = Inf, fill = Category), alpha = 0.1) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Willingness to Punish (Avg)", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +  # Only rows are faceted
  scale_x_continuous(
    limits = c(-0.41, 0.2), 
    breaks = round(seq(-0.4, 0.42, 0.2), 2)) +
  scale_y_continuous(
    breaks = 1:length(levels(punish_data$Category)),  
    labels = c("Denial", "Validity \nChallenge", "Applicability \nChallenge")) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 14, hjust = 0),  
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

punish_plot


## Figure A2: Effect of Mitigation Strategies on Economic Sanctions 

Punish <- lm(punish_econ ~ treatment, data=P_noIO)
PS <- cl(Punish, P_noIO$ResponseId)

punish_data <- data.frame(
  Estimate = coef(Punish)[-1], SE = PS[-1, 2],
  CI_low = coef(Punish)[-1] - 1.96 * PS[-1, 2],
  CI_high = coef(Punish)[-1] + 1.96 * PS[-1, 2],
  Category = factor(c("Denial", "Validity Challenge", "Applicability Challenge"),
                    levels = c("Denial", "Validity Challenge", "Applicability Challenge"))
)

punish_plot <- ggplot(punish_data, aes(y = as.numeric(as.factor(Category)), x = Estimate)) +
  geom_rect(aes(
    ymin = as.numeric(as.factor(Category)) - 0.5,
    ymax = as.numeric(as.factor(Category)) + 0.5,
    xmin = -Inf, xmax = Inf, fill = Category), alpha = 0.1) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    shape = 16, size = 0.8, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Willingness to Punish (Sanctions)", x = "", y = NULL) +
  facet_grid(rows = vars(Category), scales = "free_y", space = "free") +  # Only rows are faceted
  scale_x_continuous(
    limits = c(-0.41, 0.2), 
    breaks = round(seq(-0.4, 0.42, 0.2), 2)) +
  scale_y_continuous(
    breaks = 1:length(levels(punish_data$Category)),  
    labels = c("Denial", "Validity \nChallenge", "Applicability \nChallenge")) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 14, hjust = 0),  
    axis.text.x = element_text(size = 12),
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "none",
    strip.text.y.right = element_blank()
  )

punish_plot


## Figure A3: Effect of Mitigation Strategies by Issue Area

WTPA <- lm(punish_mean ~ treatment_expand, data=P[P$issue=="aggression",])
WTPA.SES <- cl(WTPA, P$ResponseId[P$issue=="aggression"])

WTPO <- lm(punish_mean ~ treatment_expand, data=P[P$issue=="human rights",])
WTPO.SES <- cl(WTPO, P$ResponseId[P$issue=="human rights"])

deny_data <- data.frame(
  Estimate = c(coef(WTPA)[7:9], coef(WTPO)[7:9]),
  SE = c(WTPA.SES[7:9, 2], WTPO.SES[7:9, 2]),
  CI_low = c(coef(WTPA)[7:9] - 1.96 * WTPA.SES[7:9, 2], coef(WTPO)[7:9] - 1.96 * WTPO.SES[7:9, 2]),
  CI_high = c(coef(WTPA)[7:9] + 1.96 * WTPA.SES[7:9, 2], coef(WTPO)[7:9] + 1.96 * WTPO.SES[7:9, 2]),
  Category = rep(c("Denial", "Denial + UN", "Denial + UN \n(Neg Prime)"), 2),
  Subsample = factor(rep(c("Aggression", "Human Rights"), each = 3))  # New grouping variable
)
deny_data$Category <- factor(deny_data$Category, 
                             levels = c("Denial + UN \n(Neg Prime)", "Denial + UN", "Denial"))

deny_plot <- ggplot(deny_data, aes(y = Category, x = Estimate, 
                                   color = Subsample, shape = Subsample, group = Subsample)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    size = 0.8, position = position_dodge(width = 0.4)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Denial, UN Response \n on Willingness to Punish",
       x = "", y = NULL, color = "Subsample", shape = "Subsample") +
  scale_x_continuous(
    limits = c(-0.6, 0.4), 
    breaks = round(seq(from = -0.4, to = 0.4, by = 0.1), 2)) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "bottom",  
    strip.text.y.right = element_blank()
  )

## left panel
deny_plot

applicability_data <- data.frame(
  Estimate = c(coef(WTPA)[2:4], coef(WTPO)[2:4]),
  SE = c(WTPA.SES[2:4, 2], WTPO.SES[2:4, 2]),
  CI_low = c(coef(WTPA)[2:4] - 1.96 * WTPA.SES[2:4, 2], coef(WTPO)[2:4] - 1.96 * WTPO.SES[2:4, 2]),
  CI_high = c(coef(WTPA)[2:4] + 1.96 * WTPA.SES[2:4, 2], coef(WTPO)[2:4] + 1.96 * WTPO.SES[2:4, 2]),
  Category = rep(c("Applicability", "Applicability + UN", "Applicability + UN \n(Neg Prime)"), 2),
  Subsample = factor(rep(c("Aggression", "Human Rights"), each = 3))  # New grouping variable
)
applicability_data$Category <- factor(applicability_data$Category, 
                                      levels = c("Applicability + UN \n(Neg Prime)", "Applicability + UN", "Applicability"))

applicability_plot <- ggplot(applicability_data, aes(y = Category, x = Estimate, 
                                                     color = Subsample, shape = Subsample, group = Subsample)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    size = 0.8, position = position_dodge(width = 0.4)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Applicability Challenge, UN \n Response on Willingness to Punish",
       x = "", y = NULL, color = "Subsample", shape = "Subsample") +
  scale_x_continuous(
    limits = c(-0.6, 0.4), 
    breaks = round(seq(from = -0.4, to = 0.4, by = 0.1), 2)) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "bottom",  
    strip.text.y.right = element_blank()
  )

## middle panel
applicability_plot

validity_data <- data.frame(
  Estimate = c(coef(WTPA)[10:12], coef(WTPO)[10:12]),
  SE = c(WTPA.SES[10:12, 2], WTPO.SES[10:12, 2]),
  CI_low = c(coef(WTPA)[10:12] - 1.96 * WTPA.SES[10:12, 2], coef(WTPO)[10:12] - 1.96 * WTPO.SES[10:12, 2]),
  CI_high = c(coef(WTPA)[10:12] + 1.96 * WTPA.SES[10:12, 2], coef(WTPO)[10:12] + 1.96 * WTPO.SES[10:12, 2]),
  Category = rep(c("Validity", "Validity + UN", "Validity + UN \n(Neg Prime)"), 2),
  Subsample = factor(rep(c("Aggression", "Human Rights"), each = 3))  # New grouping variable
)
validity_data$Category <- factor(validity_data$Category, 
                                 levels = c("Validity + UN \n(Neg Prime)", "Validity + UN", "Validity"))

validity_plot <- ggplot(validity_data, aes(y = Category, x = Estimate, 
                                           color = Subsample, shape = Subsample, group = Subsample)) +
  geom_pointrange(aes(
    xmin = CI_low, xmax = CI_high), 
    size = 0.8, position = position_dodge(width = 0.4)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75) +
  labs(title = "Effect of Validity Challenge, UN \n Response on Willingness to Punish",
       x = "", y = NULL, color = "Subsample", shape = "Subsample") +
  scale_x_continuous(
    limits = c(-0.6, 0.4), 
    breaks = round(seq(from = -0.4, to = 0.4, by = 0.1), 2)) +
  theme_minimal() +
  theme(
    strip.text.y = element_text(size = 14, face = "bold"), 
    axis.text.y = element_text(size = 13, hjust = 0, margin = margin(r = -20)),  
    axis.text.x = element_text(size = 12),
    plot.title = element_text(size = 15.5, face = "bold", hjust = 0.5, margin = margin(b = 15)), 
    panel.grid.major.x = element_line(color = "gray90"),
    panel.grid.major.y = element_blank(),
    panel.spacing.x = unit(2, "lines"), 
    legend.position = "bottom",  
    strip.text.y.right = element_blank()
  )

## right panel
validity_plot



## Figure A4: Legitimacy Het FX of UN interventions

WTP <- lm(punish_PC ~ treatment_expand, data=P)
BS <- cl(WTP, P$ResponseId)

WTP_NP <- lm(punish_PC ~ treatment_expand, data=P[P$IOprime=="negative",])
LPS <- cl(WTP_NP, P$ResponseId[P$IOprime=="negative"])

WTP_LC <- lm(punish_PC ~ treatment_expand, data=P[P$UN_confidence_score==0,])
LowS <- cl(WTP_LC, P$ResponseId[P$UN_confidence_score==0])

WTP_HC <- lm(punish_PC ~ treatment_expand, data=P[P$UN_confidence_score==1,])
HiS <- cl(WTP_HC, P$ResponseId[P$UN_confidence_score==1])

P$Resp_Norms_Score <- ifelse(P$Resp_Norms >=4, 1, 0)

WTP_LN <- lm(punish_PC ~ treatment_expand, data=P[P$Resp_Norms_Score==0,])
UpholdS <- cl(WTP_LN, P$ResponseId[P$Resp_Norms_Score==0])

WTP_HN <- lm(punish_PC ~ treatment_expand, data=P[P$Resp_Norms_Score==1,])
UpholdDS <- cl(WTP_HN, P$ResponseId[P$Resp_Norms_Score==1])


## Denial plot (left panel)
par(mfrow=c(1,1))
plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.4,0.3), ylim=c(6.9,8),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Deny + UN Response")
abline(v=0, lty=2)
axis(side=1, at = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2))

points(BS["treatment_expanddeny + UN",1], 8, col="black", pch=19, cex=1.35)
lines(c(BS["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*BS["treatment_expanddeny + UN",2]),
      c(8,8), lwd=1)
text(x=-0.27, y=8, "Denial + UN", cex=0.9)

points(LPS["treatment_expanddeny + UN, neg prime",1], 7.8, col="black", pch=19, cex=1.35)
lines(c(LPS["treatment_expanddeny + UN, neg prime",1] + c(-1,1)*qnorm(0.975)*LPS["treatment_expanddeny + UN, neg prime",2]),
      c(7.8,7.8), lwd=1)
text(x=-0.3, y=7.83, "Denial + UN", cex=0.9)
text(x=-0.3, y=7.77, "(Negative Prime)", cex=0.9)

points(LowS["treatment_expanddeny + UN",1], 7.6, col="black", pch=19, cex=1.35)
lines(c(LowS["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*LowS["treatment_expanddeny + UN",2]),
      c(7.6,7.6), lwd=1)
text(x=-0.3, y=7.63, "Denial + UN", cex=0.9)
text(x=-0.3, y=7.57, "(Low Conf)", cex=0.9)

points(HiS["treatment_expanddeny + UN",1], 7.4, col="black", pch=19, cex=1.35)
lines(c(HiS["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*HiS["treatment_expanddeny + UN",2]),
      c(7.4,7.4), lwd=1)
text(x=-0.26, y=7.43, "Denial + UN", cex=0.9)
text(x=-0.26, y=7.37, "(Hi Conf)", cex=0.9)

points(UpholdDS["treatment_expanddeny + UN",1], 7.2, col="black", pch=19, cex=1.35)
lines(c(UpholdDS["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*UpholdDS["treatment_expanddeny + UN",2]),
      c(7.2,7.2), lwd=1)
text(x=-0.28, y=7.23, "Denial + UN", cex=0.9)
text(x=-0.28, y=7.17, "(Low Responsibility)", cex=0.9)

points(UpholdS["treatment_expanddeny + UN",1], 7, col="black", pch=19, cex=1.35)
lines(c(UpholdS["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*UpholdS["treatment_expanddeny + UN",2]),
      c(7,7), lwd=1)
text(x=-0.28, y=7.03, "Denial + UN", cex=0.9)
text(x=-0.28, y=6.97, "(Hi Responsibility)", cex=0.9)

## Applicatory plot (middle panel)
plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.75,0.3), ylim=c(6.9,8),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Applicatory Challenge + UN")
abline(v=0, lty=2)
axis(side=1, at = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2))

points(BS["treatment_expandapplicatory + UN",1], 8, col="black", pch=19, cex=1.35)
lines(c(BS["treatment_expandapplicatory + UN",1] + c(-1,1)*qnorm(0.975)*BS["treatment_expandapplicatory + UN",2]),
      c(8,8), lwd=1)
text(x=-0.32, y=8, "Applicatory + UN", cex=0.9)

points(LPS["treatment_expandapplicatory + UN, neg prime",1], 7.8, col="black", pch=19, cex=1.35)
lines(c(LPS["treatment_expandapplicatory + UN, neg prime",1] + c(-1,1)*qnorm(0.975)*LPS["treatment_expandapplicatory + UN, neg prime",2]),
      c(7.8,7.8), lwd=1)
text(x=-0.59, y=7.83, "Applicatory + UN", cex=0.9)
text(x=-0.59, y=7.77, "(Negative Prime)", cex=0.9)

points(LowS["treatment_expandapplicatory + UN",1], 7.6, col="black", pch=19, cex=1.35)
lines(c(LowS["treatment_expandapplicatory + UN",1] + c(-1,1)*qnorm(0.975)*LowS["treatment_expandapplicatory + UN",2]),
      c(7.6,7.6), lwd=1)
text(x=-0.4, y=7.63, "Applicatory + UN", cex=0.9)
text(x=-0.4, y=7.57, "(Low Conf)", cex=0.9)

points(HiS["treatment_expandapplicatory + UN",1], 7.4, col="black", pch=19, cex=1.35)
lines(c(HiS["treatment_expandapplicatory + UN",1] + c(-1,1)*qnorm(0.975)*HiS["treatment_expandapplicatory + UN",2]),
      c(7.4,7.4), lwd=1)
text(x=-0.34, y=7.43, "Applicatory + UN", cex=0.9)
text(x=-0.34, y=7.37, "(Hi Conf)", cex=0.9)

points(UpholdDS["treatment_expandapplicatory + UN",1], 7.2, col="black", pch=19, cex=1.35)
lines(c(UpholdDS["treatment_expandapplicatory + UN",1] + c(-1,1)*qnorm(0.975)*UpholdDS["treatment_expandapplicatory + UN",2]),
      c(7.2,7.2), lwd=1)
text(x=-0.43, y=7.23, "Applicatory + UN", cex=0.9)
text(x=-0.43, y=7.17, "(Low Responsibility)", cex=0.9)

points(UpholdS["treatment_expandapplicatory + UN",1], 7, col="black", pch=19, cex=1.35)
lines(c(UpholdS["treatment_expandapplicatory + UN",1] + c(-1,1)*qnorm(0.975)*UpholdS["treatment_expandapplicatory + UN",2]),
      c(7,7), lwd=1)
text(x=-0.4, y=7.03, "Applicatory + UN", cex=0.9)
text(x=-0.4, y=6.97, "(Hi Responsibility)", cex=0.9)

## Validity plot (right panel)
par(mfrow=c(1,1))
plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.4,0.3), ylim=c(6.9,8),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Validity Challenge + UN")
abline(v=0, lty=2)
axis(side=1, at = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2))

points(BS["treatment_expandvalidity + UN",1], 8, col="black", pch=19, cex=1.35)
lines(c(BS["treatment_expandvalidity + UN",1] + c(-1,1)*qnorm(0.975)*BS["treatment_expandvalidity + UN",2]),
      c(8,8), lwd=1)
text(x=-0.27, y=8, "Validity + UN", cex=0.9)

points(LPS["treatment_expandvalidity + UN, neg prime",1], 7.8, col="black", pch=19, cex=1.35)
lines(c(LPS["treatment_expandvalidity + UN, neg prime",1] + c(-1,1)*qnorm(0.975)*LPS["treatment_expandvalidity + UN, neg prime",2]),
      c(7.8,7.8), lwd=1)
text(x=-0.34, y=7.83, "Validity + UN", cex=0.9)
text(x=-0.34, y=7.77, "(Negative Prime)", cex=0.9)

points(LowS["treatment_expandvalidity + UN",1], 7.6, col="black", pch=19, cex=1.35)
lines(c(LowS["treatment_expandvalidity + UN",1] + c(-1,1)*qnorm(0.975)*LowS["treatment_expandvalidity + UN",2]),
      c(7.6,7.6), lwd=1)
text(x=-0.27, y=7.63, "Validity + UN", cex=0.9)
text(x=-0.27, y=7.57, "(Low Conf)", cex=0.9)

points(HiS["treatment_expandvalidity + UN",1], 7.4, col="black", pch=19, cex=1.35)
lines(c(HiS["treatment_expandvalidity + UN",1] + c(-1,1)*qnorm(0.975)*HiS["treatment_expandvalidity + UN",2]),
      c(7.4,7.4), lwd=1)
text(x=-0.32, y=7.43, "Validity + UN", cex=0.9)
text(x=-0.32, y=7.37, "(Hi Conf)", cex=0.9)

points(UpholdDS["treatment_expandvalidity + UN",1], 7.2, col="black", pch=19, cex=1.35)
lines(c(UpholdDS["treatment_expandvalidity + UN",1] + c(-1,1)*qnorm(0.975)*UpholdDS["treatment_expandvalidity + UN",2]),
      c(7.2,7.2), lwd=1)
text(x=-0.28, y=7.23, "Validity + UN", cex=0.9)
text(x=-0.28, y=7.17, "(Low Responsibility)", cex=0.9)

points(UpholdS["treatment_expandvalidity + UN",1], 7, col="black", pch=19, cex=1.35)
lines(c(UpholdS["treatment_expandvalidity + UN",1] + c(-1,1)*qnorm(0.975)*UpholdS["treatment_expandvalidity + UN",2]),
      c(7,7), lwd=1)
text(x=-0.28, y=7.03, "Validity + UN", cex=0.9)
text(x=-0.28, y=6.97, "(Hi Responsibility)", cex=0.9)

## Figure A5: Treatment Effects among Foreign Policy Elites vs. US Public

Treats <- lm(punish_PC ~ treatment_expand, data=E)
coefs <- coef(Treats)[-1]
ses <- cl(Treats, E$ResponseId)[2:7,2]

P_IO <- P[-which(P$treatment=="control" & P$IOtreatment==1),]
P_IO <- P_IO[-which(P_IO$IOprime=="negative"),]

Treats2 <- lm(punish_PC ~ treatment_expand, data=P_IO)
coefs2 <- coef(Treats2)[-1]
ses2 <- cl(Treats2, P_IO$ResponseId)[2:7,2]

plot_data_combined <- data.frame(
  Estimate = c(coefs, coefs2),
  SE = c(ses, ses2),
  CI_low = c(coefs - 1.96 * ses, coefs2 - 1.96 * ses2),
  CI_high = c(coefs + 1.96 * ses, coefs2 + 1.96 * ses2),
  Category = rep(c(
    "Applicatory", "Applicatory + UN message",
    "Denial", "Denial + UN message",
    "Validity", "Validity + UN message"), 2),
  Panel = rep(c(
    "Applicatory Challenge", "Applicatory Challenge",
    "Denial", "Denial",
    "Validity Challenge", "Validity Challenge"), 2),
  Sample = rep(c("Sample 1", "Sample 2"), each = 6)
)

# Reorder factor levels for plotting
plot_data_combined$Category <- factor(
  plot_data_combined$Category,
  levels = c(
    "Denial + UN message",
    "Denial",
    "Validity + UN message",
    "Validity",
    "Applicatory + UN message",
    "Applicatory"
  ),
  labels = c(
    "Denial +\nUN message",
    "Denial",
    "Validity \nChallenge +\nUN message",
    "Validity \nChallenge",
    "Applicatory \nChallenge +\nUN message",
    "Applicatory \nChallenge"
  )
)

plot_data_combined$Panel <- factor(
  plot_data_combined$Panel,
  levels = c("Denial", "Validity Challenge", "Applicatory Challenge")
)

plot_data_combined$Position <- as.numeric(plot_data_combined$Category) +
  ifelse(plot_data_combined$Sample == "Sample 1", 0.15, -0.15)

combined_plot <- ggplot() +
  geom_rect(
    data = plot_data_combined,
    aes(
      ymin = as.numeric(Category) - 0.5,
      ymax = as.numeric(Category) + 0.5,
      xmin = -Inf,
      xmax = Inf,
      fill = Panel
    ),
    alpha = 0.1
  ) +
  geom_pointrange(
    data = subset(plot_data_combined, Sample == "Sample 1"),
    aes(
      y = Position,
      x = Estimate,
      xmin = CI_low,
      xmax = CI_high
    ),
    shape = 16, size = 0.8
  ) +
  geom_pointrange(
    data = subset(plot_data_combined, Sample == "Sample 2"),
    aes(
      y = Position,
      x = Estimate,
      xmin = CI_low,
      xmax = CI_high
    ),
    shape = 17, size = 0.8  
  ) +
  geom_vline(xintercept = 0, color = "#780000", size = 0.75) +
  facet_grid(rows = vars(Panel), scales = "free_y", space = "free") +
  labs(
    title = "Foreign Policy Elites vs. U.S. Public",
    x = "Effect on Willingness to Punish",
    y = NULL
  ) +
  scale_x_continuous(
    limits = c(-0.8, 0.47)
  ) +
  scale_y_continuous(
    breaks = 1:length(levels(plot_data_combined$Category)),
    labels = levels(plot_data_combined$Category)
  ) +
  theme_minimal() +
  theme(
    strip.text = element_blank(),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
    axis.text.y = element_text(size = 11, hjust = 0, margin = margin(r = 10)),
    axis.title.x = element_text(size = 14, margin = margin(t = 10)),
    panel.grid.major.y = element_blank(),
    panel.grid.major.x = element_line(color = "gray90"),
    axis.line = element_blank(),
    legend.position = "none"
  )

combined_plot


##################
## Appendix A.2 ##
##################

## Read in Nov 2023 public sample
O <- read.csv("Public_Nov2023.csv")
O$treatment <- relevel(as.factor(O$treatment), ref="control")
O$treatment_expand <- relevel(as.factor(O$treatment_expand), ref="control")



## Figure A.2.1: Effect of Mitigation Strategies on 
## Certainty, Outrage (Previous Study)

O_noIO <- O[O$IOtreatment==0 & O$treat_CM==0,]

Certainty_O <- lm(certainty ~ treat_expand, data=O_noIO)
cl_CO <- cl(Certainty_O, O_noIO$ResponseID)

Outrage_O <- lm(outrage ~ treat_expand, data=O_noIO)
cl_OO <- cl(Outrage_O, O_noIO$ResponseID)

## Certainty plot (left panel)
plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.65,0.36), ylim=c(6.8,7.6),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Certainty")
abline(v=0, lty=2)
axis(side=1, at = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6))

points(cl_CO["treat_expandstate deny",1], 7.5, col="black", pch=19, cex=1.35)
lines(c(cl_CO["treat_expandstate deny",1] + c(-1,1)*qnorm(0.95)*cl_CO["treat_expandstate deny",2]),
      c(7.5,7.5), lwd=1.5)
lines(c(cl_CO["treat_expandstate deny",1] + c(-1,1)*qnorm(0.975)*cl_CO["treat_expandstate deny",2]),
      c(7.5,7.5), lwd=.75)
text(x=-0.6, y=7.5, "Deny", cex=.87)

points(cl_CO["treat_expandstate attack",1], 7, col="black", pch=19, cex=1.35)
lines(c(cl_CO["treat_expandstate attack",1] + c(-1,1)*qnorm(0.95)*cl_CO["treat_expandstate attack",2]),
      c(7,7), lwd=1.5)
lines(c(cl_CO["treat_expandstate attack",1] + c(-1,1)*qnorm(0.975)*cl_CO["treat_expandstate attack",2]),
      c(7,7), lwd=.75)
text(x=-0.13, y=7, "Attack", cex=.87)

## Outrage plot (right panel)
plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.65,0.36), ylim=c(6.8,7.6),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Outrage")
abline(v=0, lty=2)
axis(side=1, at = c(-0.6, -0.4, -0.2, 0, 0.2, 0.4, 0.6))

points(cl_OO["treat_expandstate deny",1], 7.5, col="black", pch=19, cex=1.35)
lines(c(cl_OO["treat_expandstate deny",1] + c(-1,1)*qnorm(0.95)*cl_OO["treat_expandstate deny",2]),
      c(7.5,7.5), lwd=1.5)
lines(c(cl_OO["treat_expandstate deny",1] + c(-1,1)*qnorm(0.975)*cl_OO["treat_expandstate deny",2]),
      c(7.5,7.5), lwd=.75)
text(x=-0.25, y=7.5, "Deny", cex=.87)

points(cl_OO["treat_expandstate attack",1], 7, col="black", pch=19, cex=1.35)
lines(c(cl_OO["treat_expandstate attack",1] + c(-1,1)*qnorm(0.95)*cl_OO["treat_expandstate attack",2]),
      c(7,7), lwd=1.5)
lines(c(cl_OO["treat_expandstate attack",1] + c(-1,1)*qnorm(0.975)*cl_OO["treat_expandstate attack",2]),
      c(7,7), lwd=.75)
text(x=-0.48, y=7, "Attack", cex=.87)


## Figure A.2.2: Effect of Mitigation Strategies on 
## Willingness to Punish (Previous Study)

O_noIO$pastbehavior <- relevel(as.factor(O_noIO$pastbehavior), ref="noncompliant")

Pun_PCA <- lm(punish_PC ~ treatment + regimetype + pastbehavior + age + male + pid, data=O_noIO)
cl_Pun <- cl(Pun_PCA, O_noIO$ResponseID)

plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.55,0.25), ylim=c(6,7.8),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Treatment Effects on Willingness to Punish")
abline(v=0, lty=2)
axis(side=1, at = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2))

points(cl_Pun["treatmentattack",1], 7.2, col="black", pch=19, cex=1.2)
lines(c(cl_Pun["treatmentattack",1] + c(-1,1)*qnorm(0.95)*cl_Pun["treatmentattack",2]),
      c(7.2,7.2), lwd=1.5)
lines(c(cl_Pun["treatmentattack",1] + c(-1,1)*qnorm(0.975)*cl_Pun["treatmentattack",2]),
      c(7.2,7.2), lwd=.75)
text(x=-0.36, y=7.2, "Attack", cex=0.87)

points(cl_Pun["treatmentdeny",1], 7.7, col="black", pch=19, cex=1.35)
lines(c(cl_Pun["treatmentdeny",1] + c(-1,1)*qnorm(0.95)*cl_Pun["treatmentdeny",2]),
      c(7.7,7.7), lwd=1.5, col="black")
lines(c(cl_Pun["treatmentdeny",1] + c(-1,1)*qnorm(0.975)*cl_Pun["treatmentdeny",2]),
      c(7.7,7.7), lwd=.75, col="black")
text(x=-0.37, y=7.7, "Deny", cex=0.87)

points(cl_Pun["regimetypedem",1], 6.7, col="black", pch=19, cex=1.35)
lines(c(cl_Pun["regimetypedem",1] + c(-1,1)*qnorm(0.95)*cl_Pun["regimetypedem",2]),
      c(6.7,6.7), lwd=1.5, col="black")
lines(c(cl_Pun["regimetypedem",1] + c(-1,1)*qnorm(0.975)*cl_Pun["regimetypedem",2]),
      c(6.7,6.7), lwd=0.75, col="black")
text(x=-0.22, y=6.75, "Democratic", cex=0.87, col="black")
text(x=-0.22, y=6.65, "Transgressor", cex=0.87, col="black")

points(cl_Pun["pastbehaviorcompliant",1], 6.2, col="black", pch=19, cex=1.35)
lines(c(cl_Pun["pastbehaviorcompliant",1] + c(-1,1)*qnorm(0.95)*cl_Pun["pastbehaviorcompliant",2]),
      c(6.2,6.2), lwd=1.5, col="black")
lines(c(cl_Pun["pastbehaviorcompliant",1] + c(-1,1)*qnorm(0.975)*cl_Pun["pastbehaviorcompliant",2]),
      c(6.2,6.2), lwd=0.75, col="black")
text(x=-0.39, y=6.25, "History of", cex=0.87, col="black")
text(x=-0.39, y=6.15, "Compliance", cex=0.87, col="black")



## Figure A.2.3: Effect of UN Intervention on
## Willingness to Punish (Previous Study)

## Denial (left panel)
DenyO <- lm(punish_PC ~ treatment_expand, data=O)
cl_DenyO <- cl(DenyO, O$ResponseID)

plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.6,0.5), ylim=c(6.3,7.4),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Effect of Denial, UN Response \n on Willingness to Punish")
abline(v=0, lty=2)
axis(side=1, at = c(-0.4, -0.2, 0, 0.2, 0.4))

points(cl_DenyO["treatment_expanddeny",1], 7.2, col="black", pch=19, cex=1.35)
lines(c(cl_DenyO["treatment_expanddeny",1] + c(-1,1)*qnorm(0.95)*cl_DenyO["treatment_expanddeny",2]),
      c(7.2,7.2), lwd=1.5, col="black")
lines(c(cl_DenyO["treatment_expanddeny",1] + c(-1,1)*qnorm(0.975)*cl_DenyO["treatment_expanddeny",2]),
      c(7.2,7.2), lwd=.75, col="black")
text(x=-0.41, y=7.2, "Deny", cex=0.87)

points(cl_DenyO["treatment_expanddeny + UN",1], 6.6, col="black", pch=19, cex=1.35)
lines(c(cl_DenyO["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.95)*cl_DenyO["treatment_expanddeny + UN",2]),
      c(6.6,6.6), lwd=1.5, col="black")
lines(c(cl_DenyO["treatment_expanddeny + UN",1] + c(-1,1)*qnorm(0.975)*cl_DenyO["treatment_expanddeny + UN",2]),
      c(6.6,6.6), lwd=.75, col="black")
text(x=-0.215, y=6.66, "Deny +", cex=0.87, col="black")
text(x=-0.18, y=6.54, "UN response", cex=0.87, col="black")

## Attack (right panel)
AttackO <- lm(punish_PC ~ treatment_expand, data=O)
cl_AttackO <- cl(AttackO, O$ResponseID)

plot(1:3, 1:3, cex=0, bty="n",
     xlim=c(-0.6,0.5), ylim=c(6.3,7.4),
     xaxt="n", yaxt="n", ylab="", xlab="",
     main="Effect of Attack, UN Response \n on Willingness to Punish")
abline(v=0, lty=2)
axis(side=1, at = c(-0.4, -0.2, 0, 0.2, 0.4))

points(cl_AttackO["treatment_expandattack",1], 7.2, col="black", pch=19, cex=1.35)
lines(c(cl_AttackO["treatment_expandattack",1] + c(-1,1)*qnorm(0.95)*cl_AttackO["treatment_expandattack",2]),
      c(7.2,7.2), lwd=1.5, col="black")
lines(c(cl_AttackO["treatment_expandattack",1] + c(-1,1)*qnorm(0.975)*cl_AttackO["treatment_expandattack",2]),
      c(7.2,7.2), lwd=.75, col="black")
text(x=-0.41, y=7.2, "Attack", cex=0.87)

points(cl_AttackO["treatment_expandattack + UN",1], 6.6, col="black", pch=19, cex=1.35)
lines(c(cl_AttackO["treatment_expandattack + UN",1] + c(-1,1)*qnorm(0.95)*cl_AttackO["treatment_expandattack + UN",2]),
      c(6.6,6.6), lwd=1.5, col="black")
lines(c(cl_AttackO["treatment_expandattack + UN",1] + c(-1,1)*qnorm(0.975)*cl_AttackO["treatment_expandattack + UN",2]),
      c(6.6,6.6), lwd=.75, col="black")
text(x=-0.475, y=6.66, "Attack +", cex=0.87, col="black")
text(x=-0.45, y=6.54, "UN response", cex=0.87, col="black")


## Table A.2.1: Effect of Mitigation Strategies on 
## Certainty, Outrage (Previous Study)

stargazer(Certainty_O, Outrage_O,
          se = list(cl_CO[,2], cl_OO[,2]),
          dep.var.labels = c("Certainty", "Outrage"),
          covariate.labels = c("Attack", "Deny"),
          omit = c("Constant"),
          omit.stat = c("f", "ser", "rsq"))


## Table A.2.2: Effect of Mitigation Strategies on 
## Willingness to Punish (Previous Study)

stargazer(Pun_PCA,
          se = list(cl_Pun[,2]),
          dep.var.labels = c("Willingness to punish"),
          covariate.labels = c("Attack", "Deny", "Democracy", 
                               "History of Compliance"),
          omit = c("Constant", "pid", "age", "male"),
          omit.stat = c("f", "ser", "rsq"))


## Table A.2.3: Effect of UN Intervention on
## Willingness to Punish (Previous Study)

stargazer(DenyO,
          se = list(cl_DenyO[,2]),
          dep.var.labels = c("Willingness to punish"),
          covariate.labels = c("Attack", "Attack + UN", "Deny", "Deny + UN"),
          omit = c("Constant", "CM"),
          omit.stat = c("f", "ser", "rsq"))



